.ipynb .pdf

repository open issue

Binder Colab Live Code

STORMS VISUALIZATION

data = load_dac_hindcast().sel(time=slice('2018-01-29','2018-02-02'))\
    .to_dataframe().reset_index()
data.time = data.time.astype(str)
px.set_mapbox_access_token(
    'pk.eyJ1IjoianRhdXNpYSIsImEiOiJja2g0bzdxYm8wMHJrMnJsc2FkY3lqd3BxIn0.FxXNhnzy6TqYMUj4U4gjBg'
)
fig = px.scatter_mapbox(
    data,lon='longitude',lat='latitude',color='dac',animation_frame='time',zoom=2.4,
    color_continuous_scale=['purple','darkblue','blue','red','yellow','white']
)
fig.update_layout(mapbox_style="dark")
fig.update_layout(margin={"r":0,"t":0,"l":0,"b":0})
fig.show()

 loading the DAC hindcast data... 

# basic
import os, glob, sys

# arrays
import random
import numpy as np
import pandas as pd
import xarray as xr
from geopy.distance import geodesic
import great_circle_calculator.great_circle_calculator as gcc

# plotting
import matplotlib
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
import matplotlib.animation as animation
import cartopy.crs as ccrs
import cartopy.feature as fcart

# time
from datetime import datetime, timedelta
# widgets
import ipywidgets as widgets
from ipywidgets import interact, interactive, fixed, interact_manual

# append sscode to path
sys.path.insert(0, os.path.join(os.path.abspath(''), '..'))

# data viewer
from sscode.data import load_moana_hindcast, load_codec_hindcast, load_dac_hindcast, \
                        join_load_uhslc_tgs, load_geocean_tgs, load_cfsr, load_moana_hindcast_ss
from sscode.plotting.utils import plot_ccrs_nz

# load the cyclons
tcs = xr.open_dataset('../data/cyclon_tracs/IBTrACS_NZ.nc')
tcs['name'] = tcs.name.astype(str)
tcs_sel = []
for ics, name in enumerate(tcs.name.values):
    if 'NOT_NAMED' not in name:
        tcs_sel.append(ics)
tcs = tcs.isel(storm=tcs_sel)
tcs_time = []
for itc in range(len(tcs.name.values)):
    time = str(tcs.isel(storm=itc).time.dropna('date_time').values[-1])
    if int(time[:4])>1980 and int(time[:4])<2018:
        tcs_time.append(time)
tcs_time = pd.to_datetime(tcs_time).round('4D').unique()

# load the era5
slp = load_cfsr(time='2017',load_winds=(False,None))[0]

 loading the sea-level-pressure fields... 

 cropping the data to 2017 


 projected winds will not be calculated... returning the SLP... 

# load the dac
dac = load_dac_hindcast()

 loading the DAC hindcast data... 

# load the codec
codec = load_codec_hindcast()
codec_max = codec.where(codec.ss>codec.ss.quantile(dim='time',q=0.99))

 loading and plotting the CoDEC numerical data... 

# load moana
moana = load_moana_hindcast(plot=True)
moana_closest_sites = []

 loading the Moana v2 hindcast data... 

# load geotgs and uhslc
tgs = load_geocean_tgs(plot=False)
uhslctgs = join_load_uhslc_tgs(plot=True)

 loading and plotting the GeoOcean tidal guages... 


 loading and plotting the UHSLC tidal guages... 

# stations
moana_closest_sites = [1193,
 613,
 1323,
 117,
 212,
 393,
 578,
 460,
 787,
 744,
 1327,
 1319,
 1124,
 709,
 480,
 803,
 1174,
 325,
 1214,
 1193,
 613,
 1323,
 117,
 212,
 393,
 578,
 460,
 787,
 744,
 1327,
 1319,
 1124,
 709,
 480,
 803,
 1174,
 325,
 1214]

# moana.isel(site=moana_closest_sites).to_netcdf('data/SS_NZ_NunMod_MOANA.nc')

# we extract the days we have storms
storm_days = pd.to_datetime(codec_max.ss.dropna(dim='time',how='all').time.values)\
    .round('3D').unique()
storm_days = storm_days[storm_days>dac.time.values[0]] # where dac exists

def extract_time_period(storm_day):
    '''
        Extract the tie lapse
    '''
    return (str(storm_day-timedelta(days=10)), str(storm_day+timedelta(days=10)))
    
# storm day widget selection
storm_day = widgets.Dropdown(options=storm_days,
                             value=storm_days[917],
                             description='Select storm DAY: ')
storm_period_selection = interactive(extract_time_period,
                                     storm_day = storm_day)
storm_period_selection

storm_period_selection.result, storm_day.value

(('2017-03-27 00:00:00', '2017-04-16 00:00:00'),
 Timestamp('2017-04-06 00:00:00'))

# ANIMATION !!

# grids, axes and fig
fig = plt.figure(figsize=(25,22))
gs = gridspec.GridSpec(nrows=20,ncols=13)
ax1 = fig.add_subplot(gs[:9,0:3],projection=ccrs.PlateCarree(central_longitude=165))
ax2 = fig.add_subplot(gs[:9,3:6],projection=ccrs.PlateCarree(central_longitude=165))
ax3 = fig.add_subplot(gs[9:,:6],projection=ccrs.PlateCarree(central_longitude=165))
axes = [fig.add_subplot(gs[i+1,7:12]) for i in range(19)]
ax1.set_extent([140,190,-20,-60],crs=ccrs.PlateCarree())
ax2.set_extent([165,180,-30,-50],crs=ccrs.PlateCarree())
ax3.set_extent([165,180,-30,-50],crs=ccrs.PlateCarree())
gl = ax3.gridlines(crs=ccrs.PlateCarree(),draw_labels=True,
                   linewidth=2,color='gray',alpha=0.5,linestyle='--')
gl.xlabel_style = {'color':'black','size':20}
gl.ylabel_style = {'color':'black','size':20}
gl.top_labels = True
gl.left_labels = True
ax3.set_facecolor('lightblue')
fig.suptitle('--- Storm Surge viewer around New Zealand at ' \
             + str(storm_day.value)[:10] + ' (storm) ---',
             fontsize=30, y=0.92)

# initialize in range
slp_storm = slp.sel(time=slice(*storm_period_selection.result))
dac_storm = dac.sel(time=slice(*storm_period_selection.result))\
    .resample(time='1H').interpolate('nearest')
codec_storm = codec.sel(time=slice(*storm_period_selection.result))  
moana_storm = moana.sel(time=slice(*storm_period_selection.result))
tgs_storm = tgs.sel(time=slice(*storm_period_selection.result))

# points to animate
points = [ax.plot([],[],'.',c='red',markersize=20)[0] for ax in axes]

def plot_times(ax):
    for t in range(0,len(codec_storm.time.values)-20,24):
        ax.text(codec_storm.time.values[t],-3.2,
                str(codec_storm.time.values[t])[:10],
                size=18,rotation=25)

def plot_numbers(ax,sss):
    number = 0
    for lon,lat in zip(codec_storm.codec_coords_lon.values,
                       codec_storm.codec_coords_lat.values):
        if number==9 or number==2 or number==0 or number==16:
            ax.text(lon-0.6,lat-0.6,number,size=16,zorder=15,
                    backgroundcolor='white' if sss[number]<0.3 else 'red',
                    transform=ccrs.PlateCarree())
        elif number==11 or number==15:
            ax.text(lon+0.4,lat+0.6,number,size=16,zorder=15,
                    backgroundcolor='white' if sss[number]<0.3 else 'red',
                    transform=ccrs.PlateCarree())
        else:
            ax.text(lon+0.4,lat+0.3,number,size=16,zorder=15,
                    backgroundcolor='white' if sss[number]<0.3 else 'red',
                    transform=ccrs.PlateCarree())
        number += 1

# initialization function: plot the background of each frame
def init():
    slpcol = slp_storm.isel(time=0).plot.contourf(
        transform=ccrs.PlateCarree(),ax=ax1,
        cmap='bwr',levels=30,add_colorbar=False,
        vmin=101300-3000,vmax=101300+3000)
    daccol = dac_storm.isel(time=0).plot(
        transform=ccrs.PlateCarree(),ax=ax2,
        cmap='jet',vmin=-0.5,vmax=0.6,
        add_colorbar=False)
    codeccol = ax3.scatter(
        codec_storm.codec_coords_lon.values, 
        codec_storm.codec_coords_lat.values,
        s=200,c=codec_storm.isel(time=0).ss.values, 
        cmap='hot',zorder=16,label='CoDEC nodes',
        vmin=-0.2,vmax=0.6,marker='s',
        edgecolor='yellow',linewidth=3,
        transform=ccrs.PlateCarree())
    moanacol = ax3.scatter(moana_storm.lon.values,
                           moana_storm.lat.values,
                           s=50,c=moana_storm.isel(time=0).ss.values,
                           cmap='hot',zorder=14,label='Moana Dataset',
                           vmin=-0.2,vmax=0.6,transform=ccrs.PlateCarree())
    # try:
    #     ax3.scatter(tgs_storm.longitude.values, 
    #                 tgs_storm.latitude.values,
    #                 s=120,c=tgs_storm.isel(time=0).ss.values, 
    #                 cmap='hot',zorder=17,label='Real TGs',
    #                 vmin=-0.2,vmax=0.6,edgecolor='blue',
    #                 transform=ccrs.PlateCarree())
    # except:
    #     pass
    # set legend and numbers
    # ax3.legend(fontsize=16)
    plot_numbers(ax3,codec_storm.isel(time=0).ss.values)
    # plot codec/tgs time series
    for axi, ax in enumerate(axes):
        ax.plot(codec_storm.isel(name=axi).time.values,
                codec_storm.isel(name=axi).ss.values*2,
                c='royalblue',alpha=0.8,label='CoDEC')
        ilon = np.argmin(np.abs(dac_storm.longitude.values - \
            codec_storm.isel(name=axi).codec_coords_lon.values))
        ilat = np.argmin(np.abs(dac_storm.latitude.values - \
            codec_storm.isel(name=axi).codec_coords_lat.values))
        ax.plot(dac.isel(longitude=ilon,latitude=ilat).time.values,
                dac.isel(longitude=ilon,latitude=ilat).values*2,
                c='crimson',alpha=0.8,label='DAC')
        ax.plot(moana_storm.isel(site=moana_closest_sites[axi]).time.values,
                moana_storm.isel(site=moana_closest_sites[axi]).ss.values*2,
                c='orange',alpha=0.9,label='Moana Dataset')
        try:
            ax.plot(tgs_storm.isel(name=axi).time.values,
                    tgs_storm.isel(name=axi).ss.values*2,
                    c='seagreen',alpha=0.8,label='TGs')
        except:
            pass
        ax.text(codec_storm.time.values[-1],0.1,
                str(axi)+' -- '+str(codec_storm.name.values[axi]),
                size=20,backgroundcolor='white')
        if axi==0:
            ax.set_title('CoDEC (blue), DAC (red), Moana (orange) TGs (green) \n '+
                         'in different locations around NZ!!',
                         fontsize=22,loc='center',y=1.1,pad=25)
            # ax.legend(loc='center',fontsize=16,ncol=4)
        else:
            ax.set_title('')
        ax.set_ylim(-0.8,1.2)
        ax.set_xlim(codec_storm.time.values[0],codec_storm.time.values[-1])
        if axi==len(axes)-1:
            plot_times(ax)
        ax.axis('off')
        points[axi].set_data(moana_storm.time.values[0],
                             moana_storm.isel(site=axi,time=0).ss.values*2)
    for ax, cbar in zip([ax1,ax2,ax3],[slpcol,daccol,codeccol]):
        ax.coastlines(linewidth=2)
        ax.add_feature(fcart.LAND,zorder=10)
        ax.set_title('')
    
    return points  
    
# animation function.  This is called sequentially
def animate(step):
    slp_storm.isel(time=step).plot.contourf(
        transform=ccrs.PlateCarree(),ax=ax1,
        cmap='bwr',levels=30,add_colorbar=False,
        vmin=101300-3000,vmax=101300+3000)
    dac_storm.isel(time=step).plot(
        transform=ccrs.PlateCarree(),ax=ax2,
        cmap='jet',vmin=-0.5,vmax=0.6, 
        add_colorbar=False)
    ax3.scatter(codec_storm.codec_coords_lon.values, 
                codec_storm.codec_coords_lat.values,
                s=200,c=codec_storm.isel(time=step).ss.values, 
                cmap='hot',zorder=16,label='CoDEC nodes',
                vmin=-0.2,vmax=0.6,marker='s',
                edgecolor='yellow',linewidth=3,
                transform=ccrs.PlateCarree())
    ax3.scatter(moana_storm.lon.values,
                moana_storm.lat.values,
                s=50,c=moana_storm.isel(time=step).ss.values,
                cmap='hot',zorder=14,label='Moana Dataset',
                vmin=-0.2,vmax=0.6,transform=ccrs.PlateCarree())
    plot_numbers(ax3,codec_storm.isel(time=step).ss.values)
    for ax in [ax1,ax2,ax3]:
        ax.set_title('')
    ax2.set_title('Current time: '+str(slp_storm.time.values[step])[:16],
                  loc='center',fontsize=25,y=1.1,x=0.1,pad=25)
    for axi, ax in enumerate(axes):
        points[axi].set_data(moana_storm.time.values[step],
                             moana_storm.isel(site=axi,time=step).ss.values*2)
    
    return points

# call the animator.  blit=True means only re-draw the parts that have changed.
anim = animation.FuncAnimation(fig, animate, init_func=init,
                               frames=3, interval=1000, blit=True,
                               repeat_delay=5000)

from IPython.display import HTML
matplotlib.rcParams['animation.embed_limit'] = 2**32
HTML(anim.to_jshtml())

Once Loop Reflect